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ABSTRACT 

Context. The major challenges for a fully polarized radiative transfer driven approach to Zeeman-Doppler imaging are still the 
enormous computational requirements. In every cycle of the iterative interplay between the forward process (spectral synthesis) and 
the inverse process (derivative based optimization) the Stokes profile synthesis requires several thousand evaluations of the polarized 
radiative transfer equation for a given stellar surface model. 

Aims. To cope with these computational demands and to allow for the incorporation of a full Stokes profile synthesis into Doppler- 
and Zeeman-Doppler imaging applications as well as into large scale solar Stokes profile inversions, we present a novel fast and 
accurate synthesis method for calculating local Stokes profiles. 

Methods. Our approach is based on artificial neural network models, which we use to approximate the complex non-linear mapping 
between the most important atmospheric parameters and the corresponding Stokes profiles. A number of specialized artificial neural 
networks, are used to model the functional relation between the model atmosphere, magnetic field strength, field inclination, and field 
azimuth, on one hand and the individual components (/, Q, U, V) of the Stokes profiles, on the other hand. 

Results. We performed an extensive statistical evaluation and show that our new approach yields accurate local as well as disk- 
integrated Stokes profiles over a wide range of atmospheric conditions. The mean rms errors for the Stokes / and V profiles are well 
below 0.2 % compared to the exact numerical solution. Errors for Stokes Q and U are in the range of 1 %. Our approach does not 
only offer an accurate approximation to the LTE polarized radiative transfer it, moreover, accelerates the synthesis by a factor of more 
than 1000. 

Key words. Line: formation - Line: profiles - Polarization - Radiative transfer - Stars: magnetic fields - Sun: magnetic fields 



1. Introduction 

The successful application of Doppler- and Zeeman-Doppler 
imaging (DI & ZDI) in the last two decades has enormously 
contributed to our curr ent knowledge about stel l ar surface ac- 
tivity and magnetism ( Strassmeier & Rice 1998; Donati 1999; 



Piskunov & Kochukhov 2002; Strassmeiet 20021; iDonati et al.l 
2003; Kochukhov et al. 2004; Donati et al. 2006). Though it has 
its limitation, Zeeman-Doppler imaging is the only available 
method that allows to estimate the surface distribution of stel- 
lar magnetic fields. Doppler- and Zeeman-Doppler imaging rely 
on the complex and computational expensive interplay between 
the forward process of spectral line synthesis on the one hand 
and the inverse process of a regularized optimization or fitting 
procedure on the other hand. Based on a highly discretized sur- 
face model Doppler- and Zeeman-Doppler imaging iteratively 
seeks the best possible solution of the surface temperature and/or 
magnetic field surface distribution, which is compatible with the 
observed rotationally-modulated line profiles. 

The spectral synthesis plays a decisive role here and requires 
a good knowledge of accurate atomic data, continuous as well 
as line opacities, broadening mechanisms, model atmospheres 
etc. The usual way to cope with the synthesis in DI applications, 
which allows for the calculation of disk-integrated line profiles 
in a reasonable time frame, is the compilation of a precalculated 
database of local line profiles for large number of different at- 
mospheric models and lines-of-sight (LOS) (Ric e & Strassmeien 
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2000). But unlike DI, the radiative transfer in magnetized me- 
dia requires the solution of the coupled set of equation for all 
four Stokes parameters (Rees & Murphy 1987) and moreover, 
since the magnetic field is a vector quantity, a proper precalcu- 
lation quickly becomes prohibitive for ZDI applications in the 
general case where no simple global magnetic field organization 
is present. 

Due to these complexities and enormous demands many ap- 
plications of ZDI (but also DI) still rely on approximate methods 
or even completely avoid a radiative transfer modeling approach. 
Some use an ad-hoc line synthesis where Stokes / profiles are 
used from template stars or m odeled by simple Gaussian profiles 
(Don ati & Collier Cameronl ll997). from which the correspond- 
ing Stokes V profiles are derived by means of the weak- field ap- 
proximation (IDonati & Brownlll997l; IDonati et al.ll 19991) . In this 
respect, one is able to entirely bypass radiative transfer model- 
ing. 

In an effort to provide the necessary tools for a fully 
polarized radiative transfer driven approach to ZDI as well 
as to provide the diagnostic capabilities for the next genera- 
tion of a high-resolution spectropolarimeter (P EPSI) at the 8.4 
m Large Binocular Telescope (LBT) dStrassmeier et alj I2003L 
2008 ), we develope d the Zeeman-Doppler imaging code iMap 
dCarroil et alj I2007D . To accelerate the synthesis process, our 
ZDI code incorporates, besides the conventional numerical im- 
plementation of the polarized radiative transfer, a novel fast and 
accurate method for calculating local and disk-integrated syn- 
thetic Stokes profiles under local thermodynamic equilibrium 
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(LTE). This approach uses artificial neural networks (ANNs) 
to approximate the polarized radiative transfer calculations. 
Artificial neural networks have already proven to be a fast alter- 
nativ e for solar Stokes profile inversions and Zeem a n tomogra- 
phy jCarroll & StaudelpOOTt ISocas-Navarroll2005lrl Hites et al.1 
2007; Carroll & Kopf 2008). In the here presented approach 
ANNs are trained on the basis of the numerical solution of the 
polarized radiative transfer equation to find an approximate rep- 
resentation of the underlying complex mapping between atmo- 
spheric quantities like temperature, magnetic field, velocity, etc., 
and the corresponding Stokes profiles. A detailed assessment 
and evaluation of the ANN synthesis confirms the accuracy of 
our novel approach and moreover demonstrates that the ANN 
approach can accelerate the synthesis process by three orders of 
magnitude. 

This paper is organized as follows: Since our method 
strongly depends on the accuracy of the underlying numerical 
method used for the calculation of synthetic Stokes spectra, we 
first perform an in-depth analysis and assessment of our polar- 
ized radiative transfer code in Sect. [2] For this purpose, we made 
a detailed comparison between our code iMap and the COSSAM 
synthesis code (Stift 2000), which was als o part of the ben ch- 
mark test and inter-agreement analysis of IWade et al.l d200ll) . In 
Sect. |3] we give a short introduction into the ANN models we 
used in this work and describe the calculation of the training 
database for different atmospheric models, magnetic field con- 
figurations and Zeeman sensitive spectral lines. In Section|4]we 
evaluate the accuracy and performance of the trained ANNs rel- 
ative to the exact numerical solution. Sect. I4.2l presents a bench- 
mark test that demonstrates the exceptional speed of the new 
synthesis method compared to the conventional numerical pro- 
cedure. Finally, we conclude in Sect. [6] with a discussion and 
summary of the presented fast Stokes profile synthesis. 

2. Numerical polarized radiative transfer 

Because of the particular importance of an accurate Stokes pro- 
file synthesis for this work, we first give a brief overview of the 
basics of the polarized line formation under local thermal equi- 
librium (LTE) as well as a short introduction of the polarized ra- 
diative transfer module implemented within our ZDI code iMap. 
This is followed by a thorough benchmark t est with th e existing 
polarized radiative transfer code COSSAM (IStiftll2000l) . 

2.1. The transport equation 

The basic transport equation for the Stokes vector / in the fre- 
quency v reads, 

H-1 = -K V (I V -S V ), (1) 

dz 

where fi = cos (ff) and theta the angle between the surface nor- 
mal and the LOS. S Y is the source vector and K v the propagation 
matrix, which contains contributions from true absorption and 
anomalous dispersion. For the sake of a better readability, we 
omit the index for the frequency in the following. Assuming lo- 
cal thermodynamic equilibrium (LTE) and neglecting continuum 
polarization the source vector is given by 

S = B(T)e , (2) 

where B(T) is the Planck function at local temperature T and 
eo - (1, 0, 0, 0) r . The total absorption or propagation matrix is 
given by 

K = k c \ + k O . (3) 



Here, k c is the continuum absorption coefficient, kq the line ab- 
sorption coefficient, 1 denotes the 4x4 identity matrix, and O 
the line propagation matrix, given as 
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The inclination angle of the magnetic field vector relative to the 
LOS is denoted herein as y and the azimuthal angle of the field 
vector as y ■ The absorption profiles <p lh b, r and the anomalous 
dispersion profiles (/>' h r are given by the Voigt- and Faraday- 
Voigt function respectively (see, e.g jRees & Mu rphy 1987J) and 
are weighted by the normalized str ength of the respective atomic 
transitions (see, e.g JStenfloll 19941) . 



2.2. Physical foundations 

Our LTE code (the forward module of iMap) works on the ba- 
sis of a given model atmos phere, e.g., provide d by ATLAS9 
(lKuruczlll993l) or Phoenix ([H auschild t et al. 19991) or by semi- 
empirical models dHolweger & Mueller! Il974l) . The program 
reads atomic line parameters, which are extracted from the Vald 
database (Piskunov et al. 1995; Kupk a et al.lll999l) . The calcu- 
lations of the co ntinuous opacities rest on a code written by 
IWittmannl (I 19741) which accounts for free-free and bound-free 
transitions of H~, Hi, Hei, He - , Hj, H£ as well as for the 
metals C, Na, Mg. Scattering coefficients are computed for 
Rayleigh scattering by neutral hydrogen, molecular hydrogen, 
and neutral helium, and for Thomson scattering by free elec- 
trons. Radiation damping and Stark broadening are taken into 
account a s well as collisional bro adening (van der Waals) ac- 
cording to Anstee & O'Mara ( 1995). Partition functions are cal- 
culated according to [Bolton (1970) and lAller &~E verett (1972). 



Abundances, ionizati on energies, and p artition functions are sup- 
plied for 83 elements ( Wittm annl 1 9741) . The Zeeman patterns are 
calculated under L-S coupling. Voigt- and Faraday-Voigt func- 
tions are evaluated using the rational approximation method of 
iHumlicekl dl982l) . 

The numerical integration of the polarized radiative trans- 
fer is provided by a diagonal element lambda operator (DELO) 
method (Rees et al. 1989) and accounts for magneto-optical ef- 
fects. The method is implemented in its linear (original) form as 
well as in the quadratic versio n, which uses a parabo lic approxi- 
mation of the source function (ISocas-Navarro et alj |2000). 

For the disk-integration and calculation of stellar (flux) spec- 
tra, the underlying surface model in iMap is parameterized on a 
variable equal-area or equal-degree partition with a minimum el- 
ement (pixel) size of 1 x 1 degree. For each surface element, a 
local Stokes vector is calculated with respect to its position and 
underlying atmospheric parameters (abundance, Doppler veloc- 
ity, bulk velocity, temperature and pressure structures, magnetic 
field, micro- and macro turbulence). All atmospheric parame- 
ters are also allowed to vary along their local vertical direction, 
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Fig. 1. Total continuous absorption coefficient over the wavelengths (top) calculated with iMap (solid black line) and COSSAM 
(dashed red line). On the left side we see the calculation for the cold temperature regime and, on the right, for the hot temperature 
regime. In the bottom panel we see the absorption for the major contributors. Both codes show a very good agreement with rms 
error smaller than 1 %. 



which facilitates the complete description of a 3-dimensional at- 
mosphere. Limb darkening is fully accounted for by adjusting 
the depth stratification of the model atmosphere (temperature 
and pressure) for each surface element with respect to the LOS, 
i.e., recalculation of the optical depth scale according toyu (angle 
between the local normal vector and the LOS). Field structures 
on the surface (in temperature, abundance and magnetic fields) 
can either be described by setting individual surface elements or 
using spherical harmonics. 



2.3. Interagreement between iMap and COSSAM 



For th is benchmark test we follow the framework of IWade et al.l 
(2001), who performed a detailed analysis of three dif- 
feren t polarized r adiative transfer codes: CO SSAM (IStiftl 
2000) ; INVERSI O dPiskunov & Kochukhovll2002l) : ZEEMAN2 
(]Landstreetll988l) . Our goal in this section is to assess the agree- 
ment between iMap and COSSAM and by this indirectly also the 
agreement betwe en iMap and the o ther synthesis codes of the 
benchmark test of IWade et ail (120011) . 
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Fig. 2. Radiative, van-der- Waals, and quadratic Stark ef fect cal- 
culated with both codes under the asun atmosphere of lKuruca 
(1992) (solid black iMap, dashed red COSSAM). For all three 
contributions to the damping factor the two codes exhibit a very 
good agreement. 



2.3.1. Continuous opacities 

We first compare the continuous opacities and its variations over 
wavelength for different temperature regimes. We therefore cal- 
culated the continuous opacities for two different temperatures in 
the range of log(/l) = 3.5 to lo g(/l) = 4.3. W e used the Kurucz so- 
lar model atmosphere asun JKuruczlll992l) to calculate the total 
absorption coefficient at 4017 K and 6099 K which corresponds 
to a height of \og(p x ) - -2.185 and log(p x ) = 0.560 respectively. 
In Fig[T] we have plotted the resulting total continuous absorp- 
tion coefficients as calculated from both codes, COSSAM and 
iMap. Additionally, we have plotted the individual main contrib- 
utors to the continuum absorption (bound-free and free-free for 
H~ and Hi, Rayleigh and Thompson scattering). Both synthesis 
codes show a very good agreement in the continuous absorption 
coefficients as well as for the individual contributors over the 
entire wavelength range. The relative rms error is less than 1 %. 



2.3.2. Line opacities 

As line opacities are related to the total absorption and the indi- 
vidual absorption profiles, we now turn our attention to the in- 
dividual damping factors which enters into the absorption and 
anomalous dispersion profiles, i.e., Voigt - and Faraday- Voigt 
functions. We follow the investigation of IWade et alJ (12001 ) and 
compare the radiative, van-der- Waals, and quadratic Stark damp- 
ing facto rs. These are c alculated again within the asun atmo- 
sphere of Kurucz (1992). In Fig. [2] we see the run of the three 
damping factors with log(p v ). We find a good agreement be- 
tween the two codes for all three damping contributors with a 
rms error of less then 0. 1 % for the radiative damping term, 0.6 
% for the van-der- Waals term, and less then 2 % for the quadratic 
Stark effect. 
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Fig. 3. Comparison of Stokes / profiles for the zero field case 
(solid black iMap, dashed red Cossam). Both profile calcula- 
tions show a very good agreement with an rms error of less than 
0.1 % . The differences in the profiles are hardly recognizable 
in the lower plot. The top panel shows the differences AI/I C be- 
tween both profiles over the wavelength. 



2.3.3. Local line profiles - zero magnetic field 

In this section, we compare the combined effect of all physi- 
cal input parameters and compare directly the resulting line pro- 
files. For this purpose, we synthesized the local line profile of 
the Fei A 6173 A line under the asun atmosphere with zero 
magnetic field. The LOS is assumed to be normal to the surface 
(6 — 0° — » /J — cos(ff) — 1) and furthermore we have assumed 
solar iron abundance of e^> = -4.37 ( defined as \ogiiFe I n to t )• 
In Fig. [3] we show the intensity line profile calculated with both 
codes (solid black is iMap, dashed red is COSSAM). Both Stokes 
/ profiles are virtually indistinguishable within Fig. [3] The rela- 
tive rms error (relative to the continuum) between both profiles 
is less than 0.1 %, a value, which is comparable to the inter- 
agreement study of Wad e et alj d200ll) . 

2.3.4. Local line profiles - kilo Gauss magnetic field 

We now introduce a volume filling homogeneous (depth- 
independent) magnetic field of 1 kG. The inclination of the field 
(relative to the LOS) is set to 40° and the azimuth to 0°. All other 
parameters are the same as in the preceding section. Figure [4] 
again demonstrates the remarkable good agreement between the 
two codes for all Stokes parameter profiles. The rms errors are as 
follows : for the Stokes / profiles = 0.18 %; for Stokes V - 0.4 
%; for Stokes Q = 1.0 %; and for Stokes U = 0.98 %. The rms 
value s are given relative to the full amplitude (see Wad e et alj 
2001) and are well within the margins of their interagreement 
study. 

From these test calculations, we can already state that both 
synthesis codes show a very good agreement, which is compa- 
rable to the inter-agreement between the three synthesis codes 
INVERS10, ZEEMAN2, and COSSAM. 



3. The artificial neural network approach 

The basic idea of the approach is to emulate the process of 
polarized line formation by using an adaptive model, that is 
fast to evaluate, and which provides the required accuracy. 



The adaptive model we seek must have a sufficient complex- 
ity to describe the non-linear mapping of Eq. (Q]), between the 
most prominent atmospheric input parameters and the result- 
ing Stokes spectra. For this purpose, we propose a supervised 
machine learning algorithm, e.g., an artificial neural network 
(ANN) model. Successful applications of ANNs in the field of 
spectral analysis and Stokes profile inversions have been pre- 
sentedb£Carroll & Staude ( 200lll2003l) : ISocas-Navarrol (12005): 
ICarroll & Kopj d2008l) . In this section, we want to show how a 
popular type of ANN, known as multilayer perceptrons (MLPs), 
can be used for the forward modeling approach, i.e., the po- 
larized spectral line synthesis. The MLP can be regarded as a 
class of nonlinear function, which performs a multivariate map- 
ping between some input vector x and an output vector y. The 
network function represents a function composition of elemen- 
tary non-linear functions g{a). These elementary functions are 
arranged in layers whereby each of these functions in one layer 
is connected via an adaptive weight vector w to all elementary 
functions in the neighboring layers. The / - th output (compo- 
nent of the output vector y) of a three layer (of weights) MLP 
for example can be written as 



yi(x; w) = gi 



V k=0 V /=0 V j=0 



\\ 



. wfxt 



(6) 



>)) 



where x, represents the /-th component of the input vector x and 



..(i) 



w - the connecting weight from the /-th input component to the 
j-th elementary function gj in the first unit layer. The weights 
then connecting all the functions gj with the functions gu in 



w 



„(3) 



the second unit layer and w lk finally connects all functions gu 
from the second unit layer to the functions gi in the third unit 
layer. The capital letters (I,J,K) giving the numbers of elemen- 
tary functions (units) in the respective unit layer. The elemen- 
tary functions g(a), which are also called activation functions, 
are given by the following type of sigmoid function, 



(a) 



1 



1 + exp(-a) 



(7) 



The network function y(x) will thus process a given input vector 
x by propagating this vector (by multiplication with the individ- 
ual weight values and subsequent summation and evaluation of 
the different activation functions) through each layer of the net- 
work. 

The particular function that will be implemented by the MLP 
is determined by the overall structure of the network and the 
individual adaptive weight values. The process of determining 
these weight values for the network function (MLP) is called 
(supervised) training and is formulated as a non-linear optimiza- 
tion problem. This training is performed on the basis of a repre- 
sentative dataset, which includes the input to target relations of 
the underlying problem. This process is similar to a non-linear 
regression for a given data set, but as the underlying model (i.e., 
MLP) is much more general, the regression function is not re- 
stricted to a specific predetermined (or anticipated) model. In 
fact, it can be shown that a MLP provides a general fram ewor k 
for approximating arbitrary non-linear functions (Bishop 1995). 
As for all statistical learning methods particular attention must 
be given to the dataset used for the training, which should con- 
tain a representative and statistical relevant sample of the under- 
lying function we wish to approximate. In our case the input- 
to-target relation is dictated by our synthesis problem and there- 
fore the input vectors are given by the underlying atmospheric 
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Fig. 4. Comparison of Stokes / profiles (top, left), Stokes V (top, right), Stokes Q (bottom,left) and Stokes U (bottom, right) for the 
1 kG case. Calculations made with iMap are again in solid black lines and those made with COSSAM in dashed red lines. Difference 
plots are shown again on top of each profile plot. 



parameters (magnetic field, temperature etc.) and the target vec- 
tors, which reslult from the polarized radiative transfer, by the 
corresponding profiles of the Stokes spectra. The training pro- 
cess (optimization) then tries to adjust the free parameters (i.e. 
weights) of the network to find an approximation for the under- 
lying generator (function) of the training data. The optimization 
problem is formulated as a non-linear, least-square problem and 
is described by the following summed-square error function 



N C 



s=5EE(**™-*l 2 



(8) 



«=1 k=\ 



where N is the number of training patterns (i.e., pairs of input- 
target vectors) in the training set and C the number of output 
units (elements) of the network; yk{x' u ,w) the network function 
for the A:-th output component given the n-th input vector x" of 
the training sample. The value t n k is the k-th (target) component 
of the n-th training vector. 

The optimization of the network function is then performed 
by a non-linear least-square procedure. Once the network is 
successfully trained and has converged in terms of minimizing 
the error function Eq. (O for the training samples, the network 
weights are frozen and the MLP, which now represents the de- 
sired approximation of the underlying function, can be applied to 
new and unknown data coming from the same parameter domain 
as the training data. 



3.1. Synthesis of the Stokes database and principal 
component analysis 

To properly train the network function Eq. ©, we have to gen- 
erate an exhaustive training database of synthetic Stokes spec- 
tra. For the task of building the training database, we have cal- 
culated with our ZDI code iMap a large database of synthetic 
local Stokes spectra. To provide the network with the neces- 
sary information about the underlying mapping, we have iden- 
tified the following atmospheric quantities as input parameters 
for the MLP: the temperature and pressure structure of model 
atmospheres (described by the effective temperature T e ff ), gravi- 
tation (log g), iron abundance (ej? e ), the local bulk velocity {vbuik) 
of the plasma, microturbulence (v„,/ c ), macroturbulence (v mac ), 
magnetic field strength (B/ oc ), magnetic field inclination (y/ oc ), 
magnetic field azimuth ixioc) and the LOS angle between ob- 
server and the local normal (#/ oc ) of the surface element. The lat- 
ter parameter determines also the center-to-limb variation (limb 
darkening) of the disk-integrated spectra. The magnetic field pa- 
rameters are assumed to be depth independent. 

To confine the possible combination of the input parameters, 
but still provide a reasonable flexibility over a wide range of pa- 
rameters for the MLP spectrum synthesis, we divided the atmo- 
spheric quantities in a set of variable and static (fixed) param- 
eters. The static parameters are initially chosen and then fixed 
during the course of the investigation. The variable as well as 
static parameters and their associated parameter ranges are listed 
in Table Q] 
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The spectral line we have synthesized is the iron line Fe i A 
6173 A , which exhibit a normal Zeeman triplet and is a suit- 
able line for the Zeeman diagnostic because of its large effective 
Lande factor (g e ff = 2.5) and because this line experiences only 
minor contribution from surrounding blends. For each Stokes pa- 
rameter (I,Q,U,V) we have calculated 20,000 Stokes spectra by 
choosing random values for the individual input parameters from 
the listed parameter ranges (Table[TJ. The static parameters also 
enter into the calculation with the values given in Table [1] but 
were not modified. Using a random generator with a uniform 
distribution ensures moreover that the available parameter space 
was evenly sampled. Note that, even though the 5-dimensional 
parameter space for the variable input parameters is relatively 
sparsely sampled, MLPs and many ANNs in general are much 
less susceptible to what is described as the curse of dimension- 
ality (the exponential growth of hypervolume) than conventional 
parameterized methods (IBishop Ill995h . 

The spectral lines were calculated in a wavelength range of + 
1 .0 A , around the central wavelength, the spectral resolution for 
the synthesis was 0.01 A , which results in a total of 201 wave- 
length points. Since the intrinsic dimensionality of the Stokes 
profiles is usually much lower than the typical wavelength sam- 
pling, the individual line profiles can be described in a much 
more compact and reduced form than in the original wavelength 
domain JAsensio Ramos et al.ll2.007b . Taking advantage of the 
extensive training database, we have applied a principle compo- 
nent analysis (PCA) to decompose the individual Stokes profiles 
x n , by calculating the respective covariance matrix C x of our 
training database, 



C x = Yj (*»W) - *W)) (*»M) - x(A)) T 



(9) 



where A is the wavelength, n the number of individual spectral 
line profiles used for the analysis, and x the mean Stokes profile 
of all spectral lines. A new set of coordinate axes can then be 
determined by calculating the eigenvectors «,• and eigenvalues v, 
of the covariance matrix Eq. (O, 

C x s t = Vi Si . (10) 

In the PCA analysis, the eigenvectors s, are grouped in descend- 
ing order according tho their eigenvalues v,- which account for 
the degree of variance in the data (Stokes spectra). The Stokes 
spectra of the database can then be described by the orthonormal 
set of eigenvectors, 



Xk(A) = 2^ a kt iSi(A) 



(ID 



where ctk,i = Xk(A)si(A) is the scalar product (projection coeffi- 
cient) between the Stokes profile Xkiy) and the eigenvector si(A). 
Retaining only the first few dominant eigenvectors (the principal 
components), which account for the majority of the variance in 
the data, we find an efficient way of reducing the dimensionality 
of the Stokes spectra by using only the projection coefficients of 
the principal components. The representation of the k-th original 
spectra Xk can then be expressed with the first / eigenvectors as 



Xk{A) = 2_] a kjn s m (A) 



(12) 



m<l 



The projection coefficients a]c, m provide the reduced (decom- 
posed) representation of the original A:-th profile. With this tech- 
nique we have reduced the dimension of each spectra from 201 
wavelength points to 9 for Stokes / and V and 14 for Stokes Q 
and U. The rms reconstruction error for the PCA decomposed 
spectra is in all cases less then 0. 1 %. 
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Table 1. Variable and static parameters and their corresponding 
parameter ranges. 



3.2. Training of the MLP synthesis networks 

For the task of approximating the polarized line formation pro- 
cess, according to Eq. (Q} we divided the approximation of the 
Stokes profiles into 5 individual MLP network models. For each 
Stokes parameter (I,Q,U,V) we use a single MLP that calculates 
the normalized line profiles (relative to the continuum intensity), 
and one MLP network that calculates the continuum intensity on 
top of the atmosphere at the central wavelength of the respective 
spectral line. The continuum intensity is expressed in units of 



10 ergs cm s 



Even though we could have used 



one MLP model for the entire process of calculating the pro- 
files of the complete Stokes vector, we decided to split up the 
networks to avoid saturation effects during the training process, 
which may be caused by the differences in the resulting magni- 
tude of the individual Stokes parameter as well as by the different 
sensitivities (responses) to the atmospheric input parameter. 

The 5 MLP network structures are implemented with three 
layers of weights according to Eq. (O. The input and output 
structure of the network is dictated by the problem. In our case, 
each synthesis network has as much input unit (elements) as 
there are variable parameters (Table Q]), i.e., the effective tem- 
perature of the model atmosphere, the magnetic field strength, 
the inclination of the magnetic field, azimuth of the field and the 
LOS between the surface normal and the observer. The static pa- 
rameters were kept fixed to their values listed in Table Q] and do 
not appear explicitly in the network structure. The output struc- 
ture of the network is given by the decomposed representation of 
the Stokes profiles. Altogether we end up with two MLP synthe- 
sis networks (for Stokes / and V) that have 5 input parameters 
and 9 output parameters (the PCA decomposed line profiles), 
two MLP synthesis networks (for Stokes Q and U) with 5 input 
and 14 output parameters (the PCA decomposed line profiles for 
Q and U), and one MLP synthesis network for the continuum 
calculation, which has 5 input parameters and one output for the 
continuum intensity. 

Prior to the training process, the network inputs (atmospheric 
parameters) as well as the outputs (decomposed spectra) of the 
training database were rescaled. A process that is necessary for 
the input parameters due to their different magnitudes, and also 
for the output parameters due to the limited output range ([0,1]) 
of the sigmoid function (Eq. (0 in the output layer of the MLP 
networks, note that the input range is {-co, +co). Each individ- 
ual input parameter was rescaled such that all rescaled param- 
eters have zero mean and a unit standard deviation. The out- 
put values were rescaled into the linear slope [0.1,0.9] of the 
sigmoid function. The 5 synthesis networks were then trained 
on the basis of the synthetic database (20,000 Stokes profiles 
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Fig. 5. The rms error distribution between the true numerical- and the MLP-calculated Stokes profiles. The validation sample 
consists of 13,000 profiles for each Stokes parameter. Due to the asymmetry of the distribution, we also list the median besides the 
mean values. The mean error for Stokes / (top left) is 0.11 %, for Stokes V (top right) 0.17 %, and 1% for Stokes Q (bottom left) 
and Stokes U (bottom right). Note that the median values for Stokes Q and Stokes U are significantly lower than the mean values. 



for each Stokes parameter) to obtain an approximate model for 
the non-linear mapping between the atmospheric input param- 
eters and the Stokes spectra. The effective complexity of the 
network and overfitting is controlled during the training pro- 
cess by an early stopping technique (Bishop 1995). Overfitting 
is here indicated by the divergence of the network error be- 
tween the training dataset and an independent validation dataset. 
For the training process and the optimization of the MLP net- 
works we have used a c onjugate gradient algorithm (see details 
ICarroll & Staudd ( 120 01)). After experimentally testing with dif- 
ferent network topologies, in terms of numbers of hidden layers 
and units therein, we found a sufficiently good convergence and 
reduction of the squared error for the following MLPs with three 
weight layers (two hidden unit layers) : Stokes / and V 5-50-35- 
9, Stokes Q and U 5-55-55-14. 

After a successful convergence of the training process, the 
network weights of the individual MLPs are frozen and the net- 
works now (in application mode) represent the desired approxi- 
mation of the spectrum synthesis function for the iron line Fe i A 
6173 A. 



4. Results 

To evaluate the performance of the trained MLPs, we have cre- 
ated another independent validation database of 1 3,000 synthetic 
local Stokes profiles for each Stokes parameter (52,000 in total) 



with our iMap code. The combinations of the input parameters 
are again determined by a random generator from the parame- 
ter ranges listed in Table Q] The calculated Stokes spectra are 
then decomposed by the PCA technique as described above. To 
ensure a valid decomposition, the new validation spectra are de- 
composed on the basis of the training database. Also, for the 
final input scaling process ,the mean and standard deviations 
from the training database are used. To assess the accuracy of 
the MLP calculation we compare the results of the MLP net- 
work calculation with the real Stokes profiles numerically cal- 
culated with iMap using the DELO integration method. Prior to 
the comparison, the calculated (scaled) PCA profiles from the 
MLP were rescaled and converted (reconstructed) to the wave- 
length domain. 



4. 1 . Local Stokes profiles 

The impressive results for the 52,000 local test Stokes spectra 
of the validation set are illustrated in the histograms in Fig. [5] 
The relative rms error is again calculated relative to the contin- 
uum intensity for Stokes / and to the full amplitude for Stokes 
V,Q, and U. The mean rms error for the Stokes / profile calcu- 
lation with the MLP is as low as 0.1 1 %; for Stokes V, 0.17 %; 
and for Stokes Q and U, slightly above 1 %. Note that due to 
the pronounced asymmetry of the Stokes Q and U error distri- 
bution the median of the distribution is only 0.5 % . The rea- 
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son why the synthesis of the Stokes Q and Stokes U profiles 
are systematically larger than those for Stokes / and V is caused 
mainly by two effects, first, for relatively weak fields (and thus 
small profiles) Stokes Q and Stokes U are second order effects 
dLandi deglTnnocentill 19921) which results in increased relative 
error contribution, and second, both components periodically 
drop to zero for certain azimuthal and inclination values, which 
again leads to very small resulting profiles and hence to an in- 
crease in the relative error (see also discussion below). Figure[6] 
shows the distribution of the relative error compared with the nu- 
merical calculated value of the continuum intensity with iMap, 
which also show a very good agreement. 

Even though effects like the temperature variation and mag- 
netic field variation leave their characteristic imprint in the 
Stokes spectra, such that an ANN is able to approximate this 
behavior, it should be kept in mind however that this will be 
done in a complex and non-linear fashion by the transport of po- 
larized radiation. In this way, the approximation capabilities and 
accuracy of the trained MLPs are quite remarkable. The MLPs 
have clearly been able to disentangle the effects of the different 
variable input parameter to calculate the correct Stokes spectra. 
Note that for the Stokes profiles as well as for the continuum 
synthesis, the effects of the LOS angle 8, which leads to a varia- 
tion of the temperature and pressure stratification over the opti- 
cal depth scale, and which is responsible for the limb darkening 
effect, could apparently be well distinguished by the synthesis 
networks from other competing effects like the magnetic field 
strength, inclination, and the effective temperature of the model 
atmosphere. 

To display how the errors are distributed over the main input 
parameter, we show in Fig.|7]the relative rms errors over four at- 
mospheric parameters. For the Stokes / profiles (Fig. [7] top row) 
we obtain a relative homogeneous distribution for the effective 
temperature (top row, first column) as well as for the magnetic 
field azimuth (top row, fourth column) which, of course, plays 
hardly any role for the Stokes / profile. The situation is differ- 
ent for the magnetic field inclination (top row third column), 
although the rms error is quite low throughout the field incli- 
nation regime (0° to 180°) one notice the increase toward the 
longitudinal orientation of the field (close to 0° and 180°). This 
can be explained by the quasi bias (border effects) in the uni- 
formly distributed training database where the pure longitudinal 
effect with a vanishing pi-component is underrepresented in the 
database and the more common non-longitudinal effect (ellipti- 
cal polarization) with three contributions (two sigma- and one 
pi-components) is far more common. Nevertheless, note the re- 
markable agreement through the entire range. For the Stokes / 
error over the magnetic field strength (top row, second column) 
we get a somewhat counterintuitive result, which shows an in- 
crease of the error with field strength. One would assume a de- 
crease of the error since for increasing field strength the individ- 
ual Zeeman components are increasingly separated. For strong 
enough fields, such that the separation of the individual Zeeman 
components is well beyond the spectral line width, one would as- 
sume that the linear dependence of the Zeeman separation upon 
the field strength should be recognizable for the MLP and there- 
fore the functional behavior should be even easier to learn for the 
MLP. This is in principle true but the PCA decomposition of the 
spectra offsets this effect and acts in the opposite direction. The 
decomposition of a completely Zeeman-splitted Stokes / profile 
on the basis of the current dataset requires a stronger represen- 
tation of the higher-order (i.e., less significant) eigenprofiles to 
describe the entire profile correctly. This means that the projec- 
tion coefficients corresponding to these higher-order eigenpro- 
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Fig. 6. The distribution of the relative discrepancy between the 
numerically-calculated continuum intensity and the MLP calcu- 
lation. The mean rms error is 0.126 % and the median 0.061 % 
which show that the MLP is able to approximate the calculation 
of the continuum intensity with high accuracy. 



files have a significant value. Despite the simpler functional de- 
pendence in the Zeeman saturated regime, the PCA decomposed 
Stokes vector gets more complicated due to the non-negligible, 
higher-order projection coefficients. This effect could be easily 
diminished by using different datasets for the weak-field and the 
strong Zeeman saturated regime. However, note that the maxi- 
mum error for strong fields is still remarkable low (lower than 
0.22 %) and we therefore refrain from splitting the database in 
this study. 

For the Stokes V profiles (Fig. [7] second row) we also get 
very satisfying results. Over the entire temperature range (sec- 
ond row, first column) we see a low error, which is around the to- 
tal mean value of 0. 17 %. The slight peak around 5500 and 6000 
K can be explained by the coincidental accumulation of low in- 
clination angles in that region. For the rms error over the field 
strength (second row, second column) we notice a slight system- 
atic increase for smaller field strength. This can be explained 
by the smaller amplitude of the Stokes V profiles in that range 
which are harder to disentangle from other effects like temper- 
ature and inclination. For higher field strengths, the Stokes V 
profiles get increasingly separated, which again affects the PCA 
representation of these profiles and thus the reconstruction error 
(see above). The overall error for the magnetic field strength is 
also very close to the total mean of 0.17 %. As the Stokes V sig- 
nals are almost independent of the azimuthal angle (except for 
small magneto-optical effects), we see no significant variation 
along the azimuthal angle (second row, fourth column). For the 
inclination (second row, third column), we see a rise in the error 
for angles around 90° which displays the cos(y) dependence of 
the Stokes V profile on the inclination angle (Eq.|5]l, which leads 
to a rapid decrease to around 90°. Although the absolute error 
is not much affected by the decrease of the signal to around 90° 
,the relative error increases for very small signals. However, note 
that the rms error is even for the very weak signals around 90° 
lower than 0.7 %. 

For Stokes Q and U, we obtain similar good results (Fig. 
Q third and fourth row). In both cases, we see a relative homo- 
geneous error distribution over the temperature (first column) as 
well as for the field strength (second column). The error distribu- 
tion over the inclination (third column) shows a relatively strong 
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Fig. 7. The distribution of the rms errors for individual atmospheric parameters. In each row, the distributions for one Stokes 
parameter are listed: (Stokes /, first row; Stokes V, second row; Stokes Q, third row; Stokes U, fourth row). In each column, 
the distributions for a particular atmospheric parameters are given: temperature (1st column), field strength (2nd column), field 
inclination (3rd column), and the field azimuth (4th column). 



increase toward 0° and 180°. This increase can be understood 
by the smallness of the resultant Stokes Q and U profiles for 
these inclination angles and the associated increase in the rela- 
tive error. But note that even for strongly inclined fields the MLP 
synthesis yields a rms error that is smaller than 6 %. The dis- 
tribution reflects here the cos(2x) and sin(2x) dependence (see 
Eq. [5} on the azimuthal angle, which leads to a rapid decrease 
of the Stokes signal at these angles and therefore again to an 
increase of the relative error. These two effects in conjunction 
(inclination and azimuth) are the main causes for the larger rel- 
ative error of the Stokes Q and U synthesis networks compared 
to the Stokes / and V networks. But nevertheless both synthe- 
sis networks (Stokes Q and L0 have learned to approximate the 
complicated relation and perform remarkably well over the en- 
tire range of parameters. 

To obtain yet another impression about the accuracy, we have 
compiled a short list of representative Stokes I , V, Q, and U 
profiles in Fig. [8] calculated under various combinations for the 
LOS angle of the atmosphere 0, the effective temperature of the 
model atmosphere T e g, magnetic field strength \B\, the magnetic 
field inclination y and the magnetic field azimuth x- Each row 
in Fig.[8]shows the set of four Stokes profiles calculated under a 



specific atmospheric conditions. The atmospheric parameters are 
specified for each row in the Stokes / plot (first column) The pro- 
files obtained by the full numerical solution (denoted by DELO) 
are shown in solid black lines, while the MLP calculations are 
shown in dashed red lines. For all conditions (set of free param- 
eters) the plots show a very good agreement between the full 
numerical solution and the MLP calculation. 



4.2. Disk integrated Stokes profiles 

To evaluate the accumulated effect of a disk-integrated spectrum 
we have also calculated the resulting profiles for two synthetic 
test star models within the iMap code. The first model con- 
sists of a potential field where the individual contributions to 
the magnetic field vector are described by spherical harmonics 
(|Altschuler & Newkirklll969l) . For the Legendre parameter, we 
have chosen l-A and ra=3, which account for a sufficient com- 
plex surface distribution of the field strength, inclination, and az- 
imuth. For the magnetic field strength, we have assumed a peak 
value of 1400 G. The surface distribution for the radial field com- 
ponent is shown in Fig. [9] 
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Fig. 8. Each row shows the set of the four Stokes profiles (I,V,Q,U), synthesized for a particular realization of the free parameters. 
The particular values of the free parameters are always specified in the Stokes / plot (first column). Profiles calculated with the full 
numerical solution (denoted by DELO), are drawn with solid black lines, while the profiles calculated with our MLP method shown 
in dashed red lines. Difference plots for each Stokes parameter are provided again on top of each profile plot. 
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Fig. 9. The radial component of the surface field for the spherical 
harmonic test star, with 1-4, m=3. The effective temperature of 
the atmospheric model is 5000 K and the peak field strength is 
± 1400 G. 
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Fig. 10. The disk-integrated Stokes spectra for the spherical har- 
monic test star. Also, for the disk-integrated Stokes spectra, 
we see very good agreement between the numerical solution 
(DELO, solid black line) and the MLP synthesis (dashed red 
line) 



For the synthesis we used a rotational velocity with a vsini of 
26 km/s. Furthermore, we again used the iron line Fe i A 6173 A 
and a Kurucz model atmosphere (log g = 4.0) with an effective 
temperature of 5000 K. The static parameters are those listed in 
TableQ] The synthetic spectra are calculated with a surface reso- 
lution of 3 x 3 degree, which results, for the visible hemisphere, 
in 3600 surface elements. For each of these surface elements the 
underlying model atmosphere is adjusted with respect to the lo- 
cal LOS angle 6 to correctly model the center-to-limb variation 
(limb darkening). 

On the basis of the stellar model, the disk-integrated spec- 
trum is calculated within iMap, one run with the full numerical 
solution and another with the MLP synthesis module. In Fig. 
[TOl the disk-integrated Stokes spectra are shown for the spher- 
ical harmonic model. These plots show very good agreement 
between the numerical solution denoted by DELO and drawn 
with solid black lines, and the MLP synthesis shown in dashed 
red lines. No accumulating errors are visible, and the errors are 
in fact equal or even smaller than the mean errors derived from 
the statistical evaluation in Sect |4.1| The rms error is 0.09 % for 
the Stokes / profile, 0.18 % for the Stokes V profile, 0.36 % 
for the Stokes Q and 0.72 % for the Stokes U profile. To pro- 
vide another model and to test the joint effect of a temperature 
and a magnetic field variation, we construct with iMap a two 
spot model star. For the quite stellar surface, we used a Kurucz 
model atmosphere (log g = 4.0) with an effective temperature 
of 5500 K and for the two spots a model that is 1000 K cooler 
(4500 K). The maximum field strength of the spots is set to 1500 
G and the field orientation is such that the two spots have a ra- 
dial field of different polarities. The vsini is 35 km/s and the 
disk integration is performed with a surface resolution of 5 x 5 
degree. To provide a smooth continuous gradient in the temper- 
ature as well as in the magnetic field within the two spots we 
have also smeared the field and temperature distribution on the 
surface with a 2-dimensional Gaussian function. Note that in this 
model the relative weighting of the local profiles is also affected 
by the temperature variation and the corresponding continuum 
intensities, which therefore requires an accurate calculation of 
the individual continuum intensities by the MLP. The model star 
is shown in Fig. QTjwhere we have color coded the surface tem- 
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Fig. 11. A two-spot test star with a bipolar spot configuration of 
± 1500 G. The effective temperature for the quiet surface is 5500 
K and 4500 K for the spots. 



perature and depicted the field lines within the spots by blue lines 
for positive and yellow lines for a negative polarity. 

The resultant Stokes profiles for the disk integration are 
shown in Fig.Q~2l where we can see again an exceptional agree- 
ment between the numerical and the neural network calculation. 
The rms errors are 0.12 % for the Stokes / profile, 0.37 % for 
the Stokes V profile, 0.44 % for the Stokes Q, and 0.43 % for the 
Stokes U profile. 



5. Improvement in speed 

One of the main motivations for our proposed approach was of 
course to find a suitable approximate method which facilitates 
a fast synthesis of Stokes profiles. Once our MLPs are trained 
and the network weights are frozen, the network structure repre- 
sents a simple form of a composite mathematical function (see 
Sect. [3). The MLP networks can then readily be implemented 
in a quickly-evaluated programming language or stored on an 
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Fig. 12. The disk-integrated Stokes spectra for the two spot test 
star. Again, an impressive agreement between the correct numer- 
ical solution (DELO, solid black line) and the MLP synthesis 
(dashed red line) is achieved. 



even faster to evaluate hardware chip (e.g., EPROM). Both ap- 
proaches (whether software or hardware) have significant poten- 
tials to allow for an efficient parallelization. For our purpose we 
implemented the five trained MLP network structures in the pro- 
gramming language C++ on a single processor desktop com- 
puter. 

To evaluate the increase in speed provided by our ANN ap- 
proach, we made the following benchmark test. Both methods, 
the numerical integration of the polarized radiative transfer equa- 
tion with DELO (which is one of the fastest numerical methods) 
and the MLP synthesis with the MLP networks, are used to cal- 
culate 10,000 times the full set of Stokes profiles (I,Q,U,V) as 
well as the local continuum intensity. The profiles are calculated 
in a wavelength range of ±2 A around the line center with a spec- 
tral resolution of 0.01 A , which gives 401 wavelength points in 
total for each Stokes parameter and each run. The spectral line 
is again the iron line Fei A 6173 A . The DELO routine inte- 
grates the spectral lines within a depth range of log(rgooo) = 1.0 
and log(T5ooo) = -5.0. The input parameters for the numerical in- 
tegration and the MLP method are chosen randomly, but prior 
to the calculation and were kept fixed for the entire benchmark 
cycle. By this we allow the numerical integration method to cal- 
culate and compile the absorption matrix Eq. © and (0, once 
and for all in advance. Moreover, it avoids the time consuming 
reading of model atmospheres and line parameters. 

The benchmark was performed on several personal comput- 
ers with varying performance characteristics (processors ranging 
from 1 .7 GHz to 3.2 GHz and RAM from 5 12 MB to 1 GB). The 
test runs for the ANN/MLP method takes between 2.1 and 4.2 
seconds while the conventional numerical calculation of the po- 
larized radiative transfer with the DELO method takes between 
2820 to 3610 seconds. This gives an improvement in speed by 
factors of 860 to 1340 ! It is this significant improvement in 
speed - by three orders of magnitude - which facilitates a com- 
plete disk-integration of a high-resolution stellar model within a 
few seconds. With this accelerated Stokes profiles synthesis, the 
complete cycle of forward and inverse calculation for ZDI appli- 
cations now become possible within a reasonable time frame. 



6. Discussion and conclusions 

In this paper we have introduced, described and evaluated a 
novel LTE Stokes spectrum synthesis method based on ANNs. 
Artificial neural networks are used to approximate the LTE po- 
larized radiative transfer by modeling the functional relation- 
ship between the most important atmospheric parameters and 
the corresponding local Stokes spectra. The ANN approach with 
MLPs is based on an extensive training database of calculated 
synthetic Stokes profiles. Because of this crucial dependence 
for the network training we first assessed the accuracy of our 
polarized radiative transfer module integrated in our ZDI code 
iMap by a detailed comparison with the existing synthesis code 
COSSAM in Sect. [3] Both codes exhibit congruent results for sev- 
eral composite values like the continuous absorption coefficient 
and broadening parameters. Moreover, both codes show a very 
good agreement in all calculated Stokes profiles. By verifying 
the agreement between iMap and COSSAM, we also indirectly 
established an inter-agreement between iMap and INVERSE10 
as well as ZEEMAN2 , as derived from the inter-agreement study 
of lWadeetal.ld200Tl) . 

In Sect. [4] we introduced the concept of our proposed method 
and described the data preparation with a PCA method as well 
as the training process for the ANNs. A detailed evaluation of 
our ANN synthesis revealed the remarkable accuracy with which 
the ANN approach is able to approximate the polarized radiative 
transfer process. Stokes profiles of the spectral line Fe i A 6173 A 
are calculated with high accuracy for different configurations of 
the atmospheric input parameters (effective temperature, mag- 
netic field strength, field inclination, field azimuth, and LOS an- 
gle). The ANNs have clearly been able to learn the complex 
non-linear mapping between the atmospheric input parameters 
and the resulting Stokes profiles. Although we have restricted 
in this work the number of variable input parameters to 5, there 
is no principle limitation to extend the number of atmospheric 
input parameters. Moreover, the limits that we have set for the 
parameter ranges (see TableQ} are also not a principle restriction 
of the ANN approach but rather driven by the constraints to cope 
with the training and evaluation database in this study. 

The ANN approach to spectral line synthesis is also not lim- 
ited to one spectral line, the ANN method can easily be ex- 
tended to a wider wavelength range to cover several neighbor- 
ing spectral lines or successive regions of interest in the wave- 
length domain. Another interesting possibility, offered by the 
ANN approach is the modeling of co-added line profiles as used 
for polarized spect ral line extraction and reconstruction algo- 
rithms like LSD dDonatiet alj|1997h or by the P CA method 
dCarroll et alj|2007t iMartinez Gonzalez et al.ll2008l) . Instead of 
synthesizing a huge number of individual Stokes profiles to 
model, the co-added (and therefore noise reduced) mean line 
profile, an ANN synthesis can be used to directly calculate this 
co-added line profile. This would offer a vast improvement in 
speed and accuracy. 

It is also important to note that since the ANN approach 
is based on the accurate numerical modeling of line profiles, 
blends can be fully accounted for with this approach. Even heav- 
ily blended line profile can be accurately synthesized as long as 
the training database properly accounts for these blends. 

One of the main reasons to approximate or even to entirely 
bypass a full radiative transfer approach in DI or ZDI appli- 
cations, is the enormous requirement for massively calculating 
local Stokes profiles during the inversion process. In our bench- 
mark tests, we could show that our ANN approach is three orders 
of magnitudes faster then the conventional numerical integration 
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scheme. This enormous acceleration and the impressive accu- Wittmann, A. 1974, Sol. Phys., 35, 11 
racy makes the ANN synthesis approach a viable and promising 
alternative for polarized radiative transfer calculations and facil- 
itates a complete radiative transfer driven approach to DI and 
ZDI as well as for large scale Stokes profile inversions. 
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manuscript. 
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